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Abstract 

Sensitivity analysis is a classical and fundamental tool to evaluate the role of a given 
parameter in a given system characteristic. Because the phase response curve is a fun- 
damental input-output characteristic of oscillators, we developed a sensitivity analysis for 
oscillator models in the space of phase response curves. The proposed tool can be applied 
OA to high-dimensional oscillator models without facing the curse of dimensionality obstacle 

associated with numerical exploration of the parameter space. Application of this tool to a 
state-of-the-art model of circadian rhythms suggests that it can be useful and instrumental 
to biological investigations. 

(N 

, , 1 Introduction 

in 

Q Circadian entrainment is a biological process at the core of most living organisms which need 

to adapt their physiological activity to the 24 hours environmental cycle associated with earth's 
rotation (e.g. variations in light or temperature condition). This process relies on the robust 
interaction between an autonomous molecular oscillator and its environment (Fig. [l]^). Ex- 
perimental observations have shown that the system is capable to exhibit oscillations with a 
period close to 24 hours in constant environmental condition (unforced system, Fig. [Tj3) and to 
lock its oscillations (in frequency and phase) to an environmental cue with a period equal to 
t***- 24 hours (periodically forced system, Fig. [TjU). This locking phenomenon is often called (circa- 

dian) entrainment [T]. Moreover, this biological process is known to be very robust, that is, it 
maintains its performance (its period and its locking) despite internal or external perturbations 
(e.g. genetic mutations, molecular noise, variability of the environmental condition, etc.). 
t-h With recent experimental advances in biology, the molecular bases of circadian rhythms has 

been increasingly unfolded in various organisms. In most eukaryotic organisms (e.g. fungus, fly, 
or mouse), the core mechanism relies on analogous interacting positive and negative feedback 
loops with several minor alterations [2]. However, even though the architecture of those bio- 
^ logical clocks is better known, the specific design and robustness mechanisms implemented in 

those architectures remain unknown [3l 0] . 

Starting with the pioneering work of Winfree [5J [6], the Phase Response Curve (PRC) has 
emerged as a fundamental input-output characteristic of oscillators. Analogously to the static 
gain of a transfer function, the PRC measures a steady-state (asymptotic) property of the 
system response induced by an impulsive input. For the static gain, the measured property is 
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Figure 1: (A) Circadian oscillators are viewed as open dynamical systems with input u and 
output y. (B) The unforced system exhibits autonomous rhythms that occur with a period 
close to 24 hours. (C) The periodically forced system adapts the organism rhythms through 
entrainment (1:1 phase- locking) with the 24- hours stimulus associated with earth's rotation. 

the integral of the response; for the PRC, the measured property is the phase shift between 
the unperturbed and perturbed responses. Because of the periodic nature of the steady-state, 
this phase shift depends on the phase at which the system receives the impulsive input. The 
PRC is thus a curve rather than a scalar. In many situations, the PRC can be determined 
experimentally and provides unique data for the model identification of the oscillator. Likewise, 
numerical methods exist to compute the PRC from a state-space model of the oscillator. Finally, 
the PRC contains the fundamental mathematical information required to reduce a n-dimensional 
state-space model to the one-dimensional (phase) center manifold of a hyperbolic limit cycle. 

In this chapter, we review (local) sensitivity tools that provide numerical and mathematical 
grounds to the robustness analysis of oscillator state-space models in connection with experi- 
mentally available observations like the PRC or the period. We then illustrate how these tools 
can be used to make physiologically relevant predictions from mathematical models of circa- 
dian rhythms. We apply our sensitivity analysis to a state-of-the-art model [7] of 16 states and 
52 parameters and exploit the results to extract the parameters and circuits that determine the 
robustness of entrainment. 

The local proposed approach is systematic and computationally tractable. It provides a 
rapid screening of all parameters, even in high-dimensional models with a large number of 
parameters. It complements nonlocal analyses often used to assess the robustness of parameters, 
such as bifurcation analysis JH] or parameter space exploration [3] [DJ . 

The chapter is organized as follows. Section [2] reviews the notion of PRCs characterizing 
the input-output behavior of an oscillator model in the neighborhood of a stable limit cycle. 
Section [3] develops the sensitivity analysis for oscillators in terms of the sensitivity of its periodic 
orbit, its PRC, and its entrainment (phase- locking) . Section [4] provides scalar robustness mea- 
sures based on this sensitivity analysis. Section [5] illustrates how those tools permit to address 
system-theoretic questions meaningful for the robustness analysis of circadian entrainment. 

2 Open oscillator models: from state-space to phase models 

In this section, we provide a short introduction to oscillators viewed as open dynamical systems, 
that is, as dynamical systems that interact with their environment [ID]. We first recall basic 
definitions about stable periodic orbits in n-dimensional state-space models (see H2] for 
details). We then introduce (finite and infinitesimal) phase response curves as fundamental 
input-output mathematical information required for the model reduction. We finally summarize 
the standard phase reduction procedure which concentrates the phase behavior information of n- 
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Figure 2: The asymptotic phase map : #(7) — > S 1 associates with each 
point x q in the basin B(j) a scalar phase &(x g ) = 9 on the unit circle S 1 such 
that lim^+oo \\<p(t, x q , 0) — (f>(t, x p , 0) || 2 = with x p = x' y (9/uj). 



dimensional state-space models into one-dimensional phase models characterized by its angular 
frequency, its PRC, and a measurement map (see [HI [H] f° r details). 

2.1 State-space models: periodic orbits and phase maps 

We consider open dynamical systems described by nonlinear (single-input and sing le-outputQ 
time-invariant state-space models 

x = f( x ) + g{x)u, xef, u 6 M, (la) 

y = h(x), yeR, (lb) 

where the vector fields f and g, and the measurement map h support all the usual smoothness 
conditions that are necessary for existence and uniqueness of solutions. We denote by </>(•, xo, u) 



the solution to the initial value problem (la) from the initial condition xo £ M n at time 0, that 
is, (f>(0,x ,u) = x . 

An oscillator is an open dynamical system whose zero-input steady-state behavior is periodic 
rather than constant. Formally, we assume that the zero- input system x = f(x) admits a locally 
hyperbolic stable periodic orbit 7 C M. n with period T (and corresponding angular frequency uj = 
2n/T). Picking an initial condition Xq on the periodic orbit 7, this latter is described by the 
(nonconstant) T-periodic trajectory (/>(•, Xq,0) =: x 7 (-), such that x 7 (-) = x 7 (- + T). The basin 
of attraction of 7 is the maximal open set from which the periodic orbit 7 attracts. (Main 
notations are illustrated on Fig. [2j) 

Since the periodic orbit 7 is a one-dimensional manifold in M ra , it is homeomorphic to the 
unit circle S . It is thus naturally parametrized in terms of a single scalar phase. The smooth 
bijective phase map : 7 — > E 1 associates with each point x p on the periodic orbit 7 its phase 
0{x p ) =: "dp on the unit circle S , such that, 

x p - x 7 (i?p/w) = 0. (2) 

This mapping is constructed such that the image of the reference point Xq is equal to (i.e. 
<9(xq) = 0) and the progression along the periodic orbit (in absence of perturbation) produces 
a constant increase in i9. The phase variable $ : M>o — > S 1 is defined along each zero-input 
trajectory (/)(■, xq,0) starting from a point xq on the periodic orbit 7, as $(t) := 0(4>(t, xo, 0)) 
for all times t > 0. The phase dynamics are thus given by •& = uj. 

For hyperbolic stable periodic orbit, the notion of phase can be extended to any point x q 
in the basin .8(7) by defining the concept of asymptotic phase. The asymptotic phase map : 



1 For presentation convenience, we consider single-input and single-output systems. All developments are easily 
generalizable to multiple-input and multiple-output systems. 
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~~ ^ ^ associates with each point x q in the basin #(7) its asymptotic phase @(x q ) =: 9 q on 
the unit circle S , such that, 

lim \m iXq ,0)-<p(t,x~<(6 q M,0)\\ 2 = 0. (3) 

t— >-+oo 

Again, this mapping is constructed such that the image of Xq is equal to and such that the 
progression along any orbit in 6(7) (in absence of perturbation) produces a constant increase 
in 9. The asymptotic phase variable 9 : M>o — > S 1 is defined along each zero-input trajectory 
(/>(•, xo,0) starting from a point xq in the basin of attraction of 7, as 8(t) := Q(4>(t,xo,0)) for 
all times t > 0. The asymptotic phase dynamics are thus given by # = w. 

The notion of asymptotic phase variable can be extended to a nonzero-input trajectory 
(f)(-,xo,u) provided that its stays in the basin of attraction of 7. In this case, the asymptotic 
phase variable is defined as 9{t) := @(<f>(t, xq, u)) for all times t > 0. Thus the variable 0(t*) at 
an instant t* > evaluates the asymptotic phase of the point (fi(t*,xo,u) such that 

lim U{t, (j)(U,xo,u),0) -4>(t, x^(9{U)/uj), 0)|| = 0. (4) 

t— >+oo 

The asymptotic phase dynamics in the case of a nonzero input is often hard to derive. 

For presentation convenience, we introduce the map x 1 : S 1 — > 7 which associates to 
each phase 6 a point <p(9 /u , Xq , 0) = x 1 {9) on the periodic orbit. This map corresponds to 
a reparametrization of the periodic solution x 7 (-). 

The 2-7r-periodic steady-state solution 5 7 (-) and the angular frequency ui can be calculated 
by solving the boundary value problem \\b\ [TB] 

(x^)'(9) - -f{x\9)) = (5a) 

x 7 (2vr) -x 7 (0) = (5b) 
^(x 7 (0),w) = (5c) 

(where the prime • denotes the derivative with respect to 9). The boundary conditions are 



given by the periodicity condition (5b) which ensures the periodicity of the map x 7 (-) and the 
phase condition (5c) which anchors a reference position x 7 (0) along the periodic orbit. The 
phase condition tp : W 1 x M>o — > R is chosen such that it defines an isolated point on the 
periodic orbit (see |16j for details). Numerical algorithms to solve this boundary value problem 
are reviewed in [T71 Appendix]. 



2.2 Phase response curves: local information about the phase map 

For many oscillators, the structure of the asymptotic phase map is very complex. This often 
makes its analytical computation impossible and even its numerical computation intractable 
(or at least very expensive, in particular for high-dimensional oscillator models). However, in 
many situations, the global knowledge of the asymptotic phase map is not required to study 
oscillator dynamics. Instead, it is sufficient to consider a local phase information also known as 
the phase response curve. 

Starting with the pioneering work of Winfree the phase response curve of an oscillator 

has proven a useful input-output tool to study oscillator dynamics. It indicates how the timing 
of inputs affects the timing (steady-state phase shift) of oscillators. Phase response curves are 
directly related to asymptotic phase maps but capture only partial (local) information about 
them. 
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Definition 2.1. The Phase Response Curve (PRC) corresponding to an impulsive input of finite 
amplitude e (i.e. u(-) := e5(-) where S(-) is the Dirac delta function) is the map q t : S 1 — > (—it, it] 
defined as 

q t {6) := AG(x^)) = lim 9(0(t, Sf(6), eS(-))) - 9(0(t, x 7 (#), 0)) . (6) 
t^o+ <« v ' <« v ' 

post-stimulus phase pre-stimulus phase 

It associates with each point on the periodic orbit (parametrized by its phase 9) the phase shift 
induced by the input. 

In many situations, the PRC can be determined experientially (in particular for circadian 
rhythms). Moreover, it can be computed numerically by simulating the nonlinear state-space 
model and comparing the asymptotic phase shift between perturbed and unperturbed trajec- 
tories. 

A mathematically more abstract — yet very useful — tool is the infinitesimal phase response 
curve. It records essentially the same information as the finite phase response curve but for 
infinitesimally small Dirac delta input (e <C 1). 

Definition 2.2. The (input) infinitesimal Phase Response Curve (iPRC) is the map q : S 1 — > R 
defined as the directional derivative 

q (9):=DQ(x^9))[g(x'(9))] (7) 

where 

De(x)[ff\:=lim @{x + e ^- @{x) . (8) 

The directional derivative can be computed as the inner product 

De(x)[g(x)] = (V x @(x),g(x)) (9) 

where V x 0(x) is the gradient of at x. The map q x : S — > W 1 : 9 h-> V x @(x' y (9)) =: q x {9) is 
known as the state infinitesimal phase response curve. 

The (state) iPRC q x {-) can be calculated by solving the boundary value problem p^3 | fTH j [T9| , 

[201 E] 

q'M + -Ux\9)) T q x {9) = Q (10a) 
to 

q x (2<K) - q x (0) = (10b) 
{q x (9), f (x< (9))) -oj = (10c) 

(where the notation A T stands for the transpose of the matrix A) . The boundary condition ( 10b ) 
imposes the periodicity of q x {-) and the normalization condition (10c) ensures a linear increase 
at rate uj of the phase variable 9 along zero-input trajectories. Numerical methods to solve this 
boundary value problem as a by-product of the periodic orbit computation are presented in [TTJ, 
Appendix] . 

Remark 1. For small values of e (i.e. e <C 1), the PRC for impulsive input of finite amplitude 
is well approximated by the iPRC, that is, q e (-) = eq(-) + C(e 2 ). 
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Figure 3: Entrainment is studied by applying weakly connected oscillator theory to the feed- 
forward interconnection between an artificial oscillator generating the input and the actual 
oscillator. 

2.3 Phase models: entrainment 

In the weak perturbation limit, that is, for small inputs 

u(t) = eu(t), e«l, \u(t)\ < 1 for all t, (11) 

any solution (j)(t, xq, u) of the oscillator model which starts in the neighborhood of the hyperbolic 
stable periodic orbit 7 stays in its neighborhood. The n-dimensional state-space model can thus 
be approximated by a one-dimensional (continuous-time) phase model [131 IS EES [20J, [21] 

6 = 00 + eq(9)u(t) 9 € §, (12a) 

y = h{9) j/El. (12b) 

The phase model is fully characterized by its angular frequency 00 > 0, its infinitesimal phase 
response map q : S 1 — > M, and its measurement map h : S 1 — > R. 

To study entrainment through weak coupling, we can apply weakly connected oscillator 
theory [TJJ Chapter 9] by considering the input u(t) as generated by an artificial oscillator 
described by the trivial phase model 9 U = oo u , y u = h u (9 u ), where we denote by oo u the input 
angular frequency and we choose the artificial oscillator output map h u such that y u (t) = 
u(t) for all times t > 0. Moreover, the network interconnection in this case is a feedforward 
interconnection from the artificial oscillator generating the input to the studied oscillator (see 
Fig.§. 

The interconnected phase dynamics are thus given by 

9 U = oo u (13a) 

6 = u + eq(9)h u (9 u ). (13b) 

Following the weakly connected oscillator theory, we decompose the angular frequencies as 
00 = Vt + A and oo u = Q. u + A u with $7 — f2 M = 0, and the phase variables as 9 = fit + ip and 
9 U = Q u t + (p u where (p and ip u are slow phase deviations from the fast oscillations Qt and Q u t. 
The phase deviation dynamics are given by 

tp u = A u (14a) 
tp = A + eq(nt + (p)h u (Q u t + (p u ). (14b) 

Assuming that A, A n , e« 1, standard averaging techniques yield 

<Pu = &u (15a) 
ip = A + er(ip-(p u ) (15b) 

where the coupling function is given by 

1 f f 

T(^ - ipu) = lim ^ / q(Qt + <p - ip u )h u (VL u t)dt. (16) 
T Jo 
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Introducing the phase difference x = <P ~ <Pu, we have 

X = A-A u + T( X )=:V(x). 



A stable equilibrium x* of (17), that is, 



x * eS 1 :V( X *) = and V( X *) < 



correspond to a stable 1:1 phase-locking behavior (or entrainment) for (13), that is, 

6(t) - 6 u (t) = x* for all times t. 



(17) 
(18) 
(19) 



In this section, we saw that a state-space oscillator model may be reduced to a phase model 
characterized by its angular frequency (or period) and its (infinitesimal) phase response curve. 
In addition, phase models are very useful to study entrainment. It is then very natural to study 
the sensitivity of oscillators with an emphasis on those characteristics. 



3 Sensitivity analysis for oscillators 

Sensitivity analysis for oscillators has been widely studied in terms of sensitivity analysis of 
periodic orbits \22\ 123"! [2^1 [25] . Because the phase response curve is an important oscillator 
characteristic, we recently proposed a sensitivity analysis of oscillator models in the space of 
phase response curves p2]. Moreover, the sensitivity analysis in the space of PRC can be 
exploited to predict the sensitivity of the entrainment. 

We summarize those developments for nonlinear time-invariant state-space models with one 
parameteij^] 

x = f(x, A) + g(x, X)u (20a) 
y = h(x,X) (20b) 

where the constant parameter A belongs to R. 



3.1 Sensitivity analysis of a periodic orbit 

The periodic orbit 7 of an oscillator model is characterized by its angular frequency u which 
measures the 'speed' of a solution along the orbit and by the 2-7r-periodic steady-state solu- 
tion x 7 (-) which describes the locus of this orbit in the state space. The sensitivity of both 
characteristics is important. 

Given a nominal parameter value Ao, the sensitivity of the angular frequency is the scalar 
Su, € M defined as 



dX 



Ao 



lim 



Xo+h 



Ao 



h 



(21) 



where the notation *| A emphasizes the parameter value A at which the model characteristic * is 
evaluated. Likewise, the sensitivity of the 27r-periodic steady-state solution is the 27r-periodic 
function Z* : S 1 — > W 1 defined as 



dx 1 
~d\ 



(•) 



Ao 



lim 

7i->0 



IAq+A 



lAo 



h 



(22) 



2 For presentation convenience, we consider systems with a one-dimensional parameter space. All developments 
are easily generalizable to systems with a q-dimensional parameter space. 
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where the explicit dependence of the 27r-periodic steady-state solution in A is given by 

^(•)|a = 0(7Ma,*2|a,O)| a . 
From ([5]), we have, taking derivatives with respect to A, 

4(0) - -A(e)Z i (6) + \v(0)S u - -b{9) = 



U) 



UJ 



where we use the following short notations 



A(.) 


:=|(^(-),Ao), 


Ipx ■ = 


6(.) 


:=^(^(.),A ), 


i>u ■= 


s(.) 


:= /(x 7 (-),A ), 





Z x (2n) -Z x (0) = 
dtp 7 



(23) 

(24a) 

(24b) 
(24c) 

(25) 
(26) 
(27) 



Remark 2. In the literature, the sensitivity of the period is often used instead of the sensitivity 
of the angular frequency. It is the scalar St £ K defined as 



St - z 



dT 



Ao 



lim 



lA 



h 



(28) 



Both sensitivity measures are equivalent up to a change of sign and a scaling factor. The 
following relationship holds 

S T /T = -S u /u. (29) 



3.2 Sensitivity analysis of a phase response curve 

Given a nominal parameter value Ao, the sensitivity of the (input) infinitesimal phase response 
curve is the 27r-periodic function Z q : S 1 — > R defined as 



z q {-) 



dq 
dX 



(•) 



Ao 



lim 

hr+0 



?(•)! 



A +/i 



9(01 



A 



From Q, we have, taking derivatives with respect to A, 



Z q {-) = (Z qx {-),g{x"<{-), Ao)) + q x (-), ~(£ 7 0, Ao)^(-) + Ao) 



dx 



dX 



(30) 



(31) 



where the 2-7r-periodic function Z qx : S 1 — > W 1 is the sensitivity of the (state) infinitesimal phase 
response curve defined as 



Z 



qx 1 



d\ 



lim 

h-s>0 



9x(-)lAo+fc _ &0)Iao 



From (10), we have, taking derivatives with respect to A, 

z' qx {e) + -A(0) T z fo w + -c[e) T q x [e) = 

H UJ UJ 

Z ffa (2vr)-Z fo (0)=0 
(Z qx (9), v(0)) + (q x (6), Zz(9)) - S^ = 



(32) 



(33a) 

(33b) 
(33c) 
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where elements of the matrix C(-) are given by 



d 2 h 
i dxjdx k 

k=l J 

OXjOA OJ OXj 



(34) 



and where the 27r-periodic function Zy : S 1 — > M. n is the sensitivity of the vector field evaluated 
along the periodic orbit defined as 



dv 
dX 



(•) 



A 



lim 



v(-)\ 



X +h V \)\\ 



h 



(35) 



Given the explicit dependence of the 27r-periodic vector field in A 

v(-) = f(x\)\ x ,X), (36) 
we have, taking derivatives with respect to A, 

Ze(-) = g(^(-),Ao)Z 5 (-) + |{(^(-),Ao). (37) 
3.3 Sensitivity analysis of the 1:1 phase-locking 

Given a nominal parameter value Ao, the sensitivity of the phase difference x* is t ne scalar 
5 V * G K defined as 



lim 



^ I Aq+/i X I Ao 
ft 



From V(x*) = 0, we have, taking derivatives of with respect to A and using (17), 



S x * 



v'(x\ )\ 



-1 



r'(x* 



Aq/ I A 



x [Sv(x\ )] 
[SA + S r ( X *\J] 



(38) 

(39) 
(40) 



where 5 A := lim ft ^. [A|A +ft - &\\ ]/ h and S r(') '■= ^ m h^o[^ (-)\ Xo+ h ~ T (-)\\ \/ h - Considering 
that uj\\ = £1 + A|a is the sum of a parameter independent term Q and a parameter dependent 
term A, we have that S u = In addition, from (16), we have, taking derivatives with respect 
to A, 

1 



Sr(- 



lim 

T-»-+oo T 



s q {nt + -)h u (n u t)dt. 



The sensitivity of the phase difference has thus two distinct contributions: 



(41) 



(42) 



where S x *\ u := — [r'(x*|A )lAo] _1 x Sw denotes the contribution of the angular frequency sen- 
sitivity and 5" x *|r := — [r / (x*|A )U ] 1 x <SHx*|a ) denotes the contribution of the coupling 
function sensitivity at x*i the latter being closely related to the iPRC. 
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3.4 Numerics of sensitivity analysis 



Numerical algorithms to solve boundary value problems (24) and (33) are reviewed in |17} 
Appendix]. We stress that existing algorithms that compute periodic orbits and iPRCs are 
easily adapted to compute their sensitivity curves, essentially at the same numerical cost. All 
numerical tests in Section [5] have been obtained with a MATL AB numerical code available from 
the authors. 

The proposed approach is systematic and computationally tractable but it only provides 
a local sensitivity analysis in the parameter space, around a nominal set of parameter values. 
It complements more global — but less tractable — tools such as bifurcation analysis or parame- 
ter space exploration. Studying the bifurcation diagram associated with a given parameter [Bj 
or using sampling methods in the full parameter space [3j |9] are classical ways to assess the 
robustness of an oscillator: the parameter range over which the oscillation exists is a (non- 
local) indicator of the sensitivity of the oscillator to the parameters. The limitation of those 
approaches is that they are univariate (only one direction of the parameter space is explored 
in a particular bifurcation diagram) and that the exploration of the parameter space rapidly 
becomes formidable as the number of parameters grows. 



4 Scalar robustness measures for oscillators 

Testing the robustness of a model against parameter variations is a basic system-theoretic ques- 
tion. In a number of situations, the very purpose of modeling is to identify those parameters that 
influence a given system property. In the literature, robustness analysis of circadian rhythms 
mostly studies the zero-input steady-state behavior such as the period or the amplitude of 
oscillations [261 El EZj and (empirical) phase-based performance measures [281 EH1 EH [31] . In 
this section, we propose scalar robustness measures to quantify the sensitivity of the angular 
frequency, the infinitesimal phase response curve, and the 1:1 phase- locking to parameters. 

Robustness measure of the angular frequency 

The angular frequency w is a positive scalar number. The sensitivity of uj with respect to the 
parameter A is thus also a scalar number S u , leading to a scalar robustness measure R^ defined 
as 

Ru ■= \S U \ (43) 
where | • | denotes the real absolute value function. 

Robustness measure of the infinitesimal phase response curve 

In contrast, the iPRC q : S 1 — > K belongs to an infinite-dimensional space Q. The sensitivity of 
q with respect to the parameter A is thus a vector S g which belongs to the tangent space T q Q 
at q. A scalar robustness measure R q is defined as 

R q ■= \\s q \\ q = ^g g (s q ,s g ) (44) 

where ||-|| denotes the norm induced by a Riemannian metric g q (•, •) at q. In this chapter, we 
use the simplest metric for signals in £2(§ 1 )B0> that is the standard inner product, 

g q (Z q ,< q )--=(Z q ,C q )= [ U®)W)de. (45) 
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In our previous work [17], we proposed further metrics which capture equivalence properties in 
the space of phase response curves. This is motivated by the fact that, in many applications, 
it is not meaningful to distinguish among PRCs that are related by a scaling factor and/or a 
phase shift. 

Robustness measure of the 1:1 phase-locking 

The stable phase difference x* 1S a scalar phase on the unit circle S . The sensitivity of x* with 
respect to the parameter A is a scalar number S x * , leading to a scalar robustness measure R x * 
defined as 

R r :=\S x *\. (46) 

Normalized robustness measures 

When analyzing a model with several parameters (A G A C M 9 ), all robustness measures R+ 
(where * stands for any characteristic of the oscillator) are g-dimensional vectors. Each element 
of those vectors represents to the scalar robustness measure corresponding to each parameter. 
The normalized robustness measure 

R* = (47) 

\\R*\\oo 

has all its components in the unit interval [0, 1]. This normalized measure allows to rank model 
parameters according to their relative ability to influence the characteristic *. 

Remark 3. Note that Rt = Rw- 

5 Application to a model of circadian rhythms 

We illustrate our sensitivity analysis on the genetic oscillator model of Leloup and Goldbeter 
(see Fig. [4]) . This model accounts for several regulatory processes identified in circadian rhythms 
of mammals. A negative autoregulatory feedback loop established by the per (period) and cry 
(cryptochrome) genes is at the heart of the circadian oscillator. The PER and CRY proteins 
form a complex PER-CRY that indirectly represses the activation of the Per and Cry genes. 
The PER-CRY complexes exert their repressive effect by binding to a complex of two proteins 
CLOCK-BMAL1. This latter, formed by the products of Clock and Bmall genes, activates 
Per and Cry transcription. In addition to this negative autoregulation, an (indirect) positive 
regulatory feedback loop is also involved. Indeed, the Bmall expression is subjected to negative 
autoregulation by CLOCK-BMAL1, through the product of the Rev-Erba gene. The complex 
PER-CRY enhances Bmall expression in an indirect manner by binding to CLOCK-BMAL1, 
and thereby reducing the transcription of the Rev-Erba gene. Finally, environmental periodic 
cycles associated with earth's rotation are mediated through light-dark cycles. Light acts on 
the system by inducing the expression of the Per gene. 

The detailed computational model of Leloup and Goldbeter possesses 16 state variables and 
52 parameters. State-space model equations and nominal parameter values are available in [71 
Supporting Text]. The effect of light is incorporated through periodic square- wave variations 
in the maximal rate of Per expression (i.e. the value of the parameter v s p goes from a constant 
low value during dark phase to a constant high value during light phase). Parameters values 
remain to be determined experimentally and have been chosen semiarbitrarily in physiological 
ranges in order to satisfy experimental observations. 
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Figure 4: The Leloup-Goldbeter model accounts for several regulatory processes identified in 
circadian rhythms of mammals. Reproduction of a figure from [TJ. 
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Each parameter of the model describes a single regulatory mechanism such as transcription 
and translation control of mRNAs, degradation of mRNAs or proteins, transport reaction, and 
phosphorylation/dephosphorylation of proteins. The analysis of single-parameter sensitivities 
reveals thus the importance of individual regulatory processes on the function of the oscillator. 

However, in order to enlighten the potential role of circuits rather than single-parameter 
properties, we grouped model parameters according to the mRNA loop to which they belong. 
Each group of parameters is associated with a different color: Per-loop in blue, Cry-loop in 
red, and Bmall -loop in green. In addition, we gathered parameters associated with interlocked 
loops in a last group represented in gray. 

In the following, we consider sensitivities to relative variations of parameters. We write 
without distinction about period sensitivities and angular frequency sensitivities due to their 
direct proportional relationship (see remarks^ and^. 

Sensitivity analysis of the period and the phase response curve 

The period and the PRC are two intrinsic characteristics of the circadian oscillator with physi- 
ological significance. We use the sensitivity analysis of the period and the PRC to measure the 
influence of regulatory processes on tuning the period and shaping the PRC. 

A two-dimensional (R u , Rq) scatter plot in which each point corresponds to a parameter of 
the model reveals the shape and strength of the relationship between both normalized robustness 
measures R^ (angular frequency or, equivalently, period) and R q (PRC). It enables to identify 
which characteristic is primarily affected by perturbations in individual parameters: parameters 
corresponding to points situated below the dashed bisector influence mostly the period; those 
above the dashed bisector influence mostly the PRC (see Fig. [5]) . 

At a coarse level of analysis, the scatter plot reveals that most parameters exhibit both low 
period and PRC sensitivities (most points are close to the origin); only few parameters display 
a medium or high sensitivity either to period or to PRC. 

At a finer level of analysis, the scatter plot reveals that the parameters associated with each 
of the three mRNA loops have distinct sensitivities: 

• the Bmall -loop parameters are associated with a high period sensitivity and a medium 
PRC sensitivity (regression line below the bisector); 

• the Per-loop parameters are associated with a medium period sensitivity and a high PRC 
sensitivity (regression line above the bisector); 

• the Cry-loop parameters are associated with a low period sensitivity and a high PRC 
sensitivity (regression line above the bisector, close to the vertical axis). 

In each feedback loop, the three more sensitive parameters represent the three same biological 
functions: the maximum rates of mRNA synthesis (v s p, v s p, and v s c), the maximum rate of 
mRNA degradation (i> m B, v m p, and v m c), and the inhibition (I) or activation (A) constants for 
the repression or enhancement of mRNA expression by BMAL1 (K~ip, Kj^p, and -KTac)- 

The small number of highly sensitive parameters is in agreement with the robust nature of 
the circadian clock and the concentration of fragilities in some specific locations of the archi- 
tecture [3]. Our analysis suggests that the transcriptional and translational control of mRNA 
(i.e. the control of both biological steps required to synthesize a protein) has to be regulated 
by specific mechanisms (not included in the model) in order to avoid failures in the clock func- 
tion. While the topology of Per- and Cry-loops are identical, the asymmetry introduced by 
the choice of parameter values leads to different sensitivity for those loops. Both loops have a 
similar high sensitivity of the PRC (while the light acts only on the maximum rate of Per mRNA 



13 







1 



Figure 5: Normalized robustness measures R w (angular frequency) and R q (iPRC) reveal the 
distinct sensitivity of three distinct genetic circuits ( Cry, Per, and Bmall ) . Each point is asso- 
ciated to a particular parameter. The three lines are regression over the parameters of the three 
gene loops. The dashed bisector indicates the positions at which both measures of robustness 
are identical. Only parameters associated with the Cry-loop exhibit low angular frequency and 
high iPRC sensitivities. The color code corresponds to different subsets of parameters associated 
to different loops (see the text for details). 
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Figure 6: Normalized sensitivity measures S x * / II •S'x* II oo ( en trainment) are due to two contribu- 
tions: S x *\ u / \\S X * (angular frequency) and S x *^/ \\S X * (coupling function). Each (thick) 
horizontal bar corresponds to a sensitivity measure with respect to a particular parameter. The 
(thin) horizontal lines indicate (in absolute value) the maximal sensitivity (among all parame- 
ters) and may be useful to compare the sensitivity of a parameter to the maximal sensitivity. 
The color code corresponds to different subsets of parameters associated to different loops (see 
the text for details). 



synthesis) but a different sensitivity of the period, the Per- loop being more sensitive than the 
Cry-loop. The high sensitivity of the period for parameters associated with the Bmall -loop 
has also being identified in [8]. However, this last prediction of the model (high sensitivity 
of the period to Bmall -loop) is not in agreement with experimental observations in [32} 133] . 
This observation may encourage the biologist and the modeler to design of new experiments to 
enlighten biological mechanisms responsible for this discrepancy between the experiment and 
the model. 



Sensitivity analysis of the entrainment 



Entrainment is an important characteristic of the circadian model. In Section 3.3, we have seen 
that the entrainment sensitivity S x * is mathematically given by the summation of two terms: a 
term S x *\ u proportional to the period sensitivity and a term S x *\r proportional to the coupling 
function sensitivity at x* ■ Those two terms correspond to two biologically distinct mechanisms 
by which the entrainment properties of the circadian clock can be regulated: a modification of 
the period or a modification of the coupling function (resulting from the modification of the 
iPRC or the input signal). 

Bar plots of S x * / \\S X * H^, S x *\ u / \\S X * H^, and S x *\y/ \\S x * in which each bar corresponds 
to a parameter allows to identify the most sensitive parameters for entrainment and to quantif}]^] 



3 The entrainment sensitivity and the contributing terms are normalized by ||S X * || (the same maximal value of 
the entrainment sensitivity) such that the summation of normalized terms is equal to the normalized entrainment 
sensitivity. 
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the respective contribution of both mechanisms in the entrainment sensitivity (see Fig. [6]) . For 
each bar plot, we sorted parameters by absolute magnitude and restricted the plot to the 14 
parameters with the highest sensitivity measure (the number 14 results from our choice to 
keep the parameters with an entrainment sensitivity greater than 0.1). Those plots allow to 
identify the parameters which play an important role in the entrainment sensitivity. We note 
that the parameter orders for S x * / ||<S' x *|| oo and S x *\ u / HS^Hoo are a l mos t identical, except for 
parameters associated with the Cry-loop. Those parameters appear in the highest ones for 

S x*\t/W S x*\\oo- 

Figure [7] (top) reveals the competitive and complementary nature of both contributions to 
entrainment sensitivity. For most parameters, both contributions have opposite signs, that is, 
points are located in the second and fourth quadrants. In addition, both mechanisms are well 
decoupled such that, when one mechanism is active, the other is almost inactive (points are 
located close to the horizontal and vertical axes). Parameters associated with Cry-loop seem 
to influence the entrainment sensitivity through a modification of the coupling function (points 
close to the vertical axis); others parameters associated with Per-loop and Bmall -loop seem to 
influence the entrainment sensitivity through a modification of the period (points close to the 
horizontal axis). 

The different mechanisms leading to entrainment sensitivity are also observed in both other 
scatter plots (see Fig. [7] bottom-left and -right). In those plots, parameters associated with 
points close to the bisector of the first and third quadrants influence the entrainment sensitivity 
through a modification of the period (bottom-left) or the coupling function (bottom- right), re- 
spectively. Again, only parameters associated with the Cry-loop seem to affect the entrainment 
through a variation of the PRC. 

Two of the parameters belonging to the Cry-loop (with high coupling function and low 
period sensitivities) have been identified by numerical simulations as important for entrainment 
properties of the model without affecting the period: Kac in [Z] and v m c in [8]. Our approach 
supports the importance of those two parameters and identifies the potential importance of a 
third one (v s c)- 

We stress that the sensitivity analysis in (TIE] is a global approach that relies on exploring the 
parameter space through numerical simulations of the model to determine the system behavior 
under constant and periodic environmental conditions while varying one parameter at a time. 
In contrast, the proposed analysis is local but systematic and computationally tractable. In the 
particular model studied here and in [8] , the predictions of the (local) sensitivity analysis match 
the predictions of the (nonlocal) analysis. 

To evaluate the nonlocal nature of our local predictions, we plot in Fig. [8] the time behavior 
of solutions for different finite (nonlocal) parameter changes. The left plots illustrate the au- 
tonomous oscillation of the isolated oscillator whereas the right plots illustrate the steady-state 
solution entrained by a periodic light input. Parameter perturbations are randomly taken in a 
range of ±10% around the nominal parameter value. Each panel corresponds to the perturba- 
tion of a different group of parameters (the black time-plot corresponds to the nominal system 
behaviors for nominal parameter values). 

A. Perturbations of three most sensitive parameters of Cry-loop (v s c, f m c> and Kxc) lead to 
small variations (mostly shortening) of the autonomous period and (not structured) large 
variations of the phase-locking. This observation is consistent with the low sensitivity of the 
period and the high sensitivity of the PRC. 

B. Perturbations of three most sensitive parameters of Bmall -loop (v s b, w m B, and K\-q) lead to 
medium variations of the autonomous period and medium variations of the phase-locking. 
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Figure 7: Normalized sensitivity measures S x * / \\S X * (entrainment), S' x *[ a ,/ \\S X * {{^ (angu- 
lar frequency), and S x *|r/ \\£>x* Hoc ( cou P nn g function) exhibit particular correlation shapes. 
The top graph represents the (S x *^/ \\S X * , >5 x *|r/ \\S X * || 00 )-plan; the bottom-left graph rep- 
resents the (5^*1^/ \\S X * , S x * / \\S X * || oo )-plan; and the bottom-right graph represents the 
(S x *\r/ ll'S'x'lloo )^x*/ ll^x* lloo)"P^ an - Each point is associated to a particular parameter. The 
color code corresponds to different subsets of parameters associated to different loops (see the 
text for details) . Those correlations support the competitive nature of both mechanisms (mod- 
ification of the period or the coupling function) leading to the entrainment sensitivity. 
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The variations of the phase-locking exhibit the same structure as variations of the period, 
suggesting that the change in period is responsible for the change of phase-locking for those 
parameters. This observation is consistent with the high sensitivity of the period and the 
medium sensitivity of the PRC. 

C. Perturbations of three most sensitive parameters of Per-loop (v s p, u m p, and Kap) exhibit 
an intermediate behavior between the situations A and B. 

D. Perturbations of parameters of interlocked loops lead to small variations of the autonomous 
period and the phase-locking, which is consistent with their low sensitivity. 

Those (nonlocal) observations are thus well predicted by the classification of parameters sug- 
gested by the (local) sensitivity analysis (see Fig. [5]). 

6 Conclusion 

This chapter proposes (local) sensitivity tools to analyze oscillator models as open dynamical 
systems. We showed that, under the weak perturbation assumption, state-space models can 
be reduced to phase models characterized by their angular frequency and their phase response 
curve. Those phase models are then useful to study the entrainment (or phase-locking) to a 
periodic input. We then introduced the sensitivity analysis for oscillators and their phase-locking 
behavior. 

The application of this approach to a detailed computational model of circadian rhythms 
provides physiologically relevant predictions. It enlightens the distinct role of different circuits in 
the robustness of entrainment and it selects 3 out of 52 parameters as parameters that strongly 
affect the phase response curve while barely affecting the period. The importance of two of 
these parameters was previously identified in the literature through simulations of the model. 

7 Lessons learnt 

Sensitivity analysis is a classical and fundamental tool to evaluate the role of a given parameter 
in a given system characteristic. Because the phase response curve is a fundamental input- 
output characteristic of oscillators, we developed a sensitivity analysis for oscillator models in 
the space of phase response curves. The proposed tool can be applied to high-dimensional 
oscillator models without facing the curse of dimensionality obstacle associated with numerical 
exploration of the parameter space. Application of this tool to a state-of-the-art model of 
circadian rhythms suggests that it can be useful and instrumental to biological investigations. 
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